ERanalysis
JMH
10/27/2020
Analysis of GPP response to N and P enrichment
Libraries
library(tidyverse)
library(gamm4)
library(mgcViz)
library(MuMIn)
library(forecast)
library(gridExtra)
library(tidymv)Load & munge data
Loads data, calculates log GPP and ER, does a ridiculous transformation of NEP, and centers stream water temp at 10C as inverse kT
met2 <- read.csv("00_metab_for_Hoodie.csv") %>%
mutate(Pdt = as.POSIXct(Pdt, format = "%m/%d/%y"),
doy = as.numeric(strftime(Pdt, format = "%j")),
Y = as.numeric(strftime(Pdt, format = "%Y")),
Yf = as.factor(as.character(Y)),
lGPP = log(GPP_median),
lER = log(-ER_median),
#centering
invKT.C = one_over_kT - 1/((10 + 273.15)*0.0000862),
streamTreat = as.factor(paste0(stream,"_",treatment))) %>%
mutate_at(vars(stream, treatment), funs(as.factor(.))) %>%
rename(invKT = one_over_kT)For predictor variables, we want to center Temp, light, Q, and GPP. We’ll do that using stream means, which we calculate below
# stream mean
metS <- met2 %>%
group_by(stream, treatment) %>%
summarise_at(vars(meanTemp, invKT.C, meanQds, LightPerDay, lGPP), funs(mean(., na.rm = TRUE))) %>%
group_by(stream) %>%
summarise_at(vars(meanTemp, invKT.C, meanQds, LightPerDay, lGPP), funs(mean(., na.rm = TRUE))) %>%
rename_at(vars(-stream), function(x) paste0(x,".StMean"))
# global mean
metG0 <- met2 %>%
group_by(stream, treatment) %>%
summarise_at(vars(meanTemp, invKT.C, meanQds, LightPerDay, lGPP), funs(mean(., na.rm = TRUE))) %>%
group_by(stream) %>%
summarise_at(vars(meanTemp, invKT.C, meanQds, LightPerDay, lGPP), funs(mean(., na.rm = TRUE))) %>%
group_by() %>%
summarise_at(vars(meanTemp, invKT.C, meanQds, LightPerDay, lGPP), funs(mean(., na.rm = TRUE)))%>%
rename_at(vars(meanTemp:lGPP), function(x) paste0(x,".GMean"))
metG <- rbind(metG0, metG0, metG0, metG0)
metG$stream <- c("st11U", "st18", "st6", "st9")Now, we log transform Q, GPP, and light then center on mean logged value (you can’t center then log).
met <- met2 %>%
full_join(metS, by = "stream") %>%
full_join(metG, by = "stream") %>%
#these are centered on the stream mean
mutate(Tanom.iktCs = invKT.C - invKT.C.StMean,
Tanom.iktCg = invKT.C - invKT.C.GMean,
Tanom.Cs = meanTemp - meanTemp.StMean,
Tanom.Cg = meanTemp - meanTemp.GMean,
LightPerDayCs = LightPerDay - LightPerDay.StMean,
LightPerDayCg = LightPerDay - LightPerDay.GMean,
meanQdsCs = meanQds - meanQds.StMean,
meanQdsCg = meanQds - meanQds.GMean,
LightPerDay.l = log(LightPerDay),
meanQds.l = log(meanQds),
LightPerDay.lcs = LightPerDay.l - log(LightPerDay.StMean),
LightPerDay.lcg = LightPerDay.l - log(LightPerDay.GMean),
meanQds.lcs = meanQds.l - log(meanQds.StMean),
meanQds.lcg = meanQds.l - log(meanQds.GMean),
lGPP.cs = lGPP - lGPP.StMean,
lGPP.cg = lGPP - lGPP.GMean,
streamTreat = paste0(stream," ", treatment),
stream = as.factor(stream))Quick visualization
ER
Focus on ER response.
Access normality: log ER - not terrible
Light
Untransformed daily integrated light
log-centered daily integrated light - stream
log-centered daily integrated light - global
Is light normal? Probably close enough.
Discharge
Untransformed mean daily Q.
log-centered mean daily Q -stream 
log-centered mean daily Q - global 
Is it normal. Log centered mean daily Q.
Temperature as 1/kT
Untransformed temperature
Temperature anom based on mean stream temp
Temperature anom based on mean stream temp - global
Normal? - looks awesome!
GPP
Is really important and there’s super interesting things going on between ER and GPP, but I’m going to leave out.
ggplot(met, aes(y = lER, x = lGPP.cg, color = stream, shape = treatment)) +geom_point() +facet_grid(stream~.)
Normal? - S18 has particularly long tails 
ggplot(met, aes(y = log(-ER_median), x = meanQds.lcg, color = Yf)) +
geom_point() +
facet_grid(. ~ stream)
Ugh - Q and temp are correlated These can’t be in the same model
ggplot(met, aes(y = meanQdsCg, x = Tanom.iktCg, color = stream)) +
geom_point() +
facet_grid(stream ~.)
met %>%
group_by(stream) %>%
summarize(currelation = cor(meanQdsCg, Tanom.iktCg))## `summarise()` ungrouping output (override with `.groups` argument)## # A tibble: 4 x 2
## stream currelation
## <fct> <dbl>
## 1 st11U 0.517
## 2 st18 0.714
## 3 st6 0.694
## 4 st9 0.487ER v. temp
ggplot(met, aes(y = lER, x = Tanom.iktCg, color = Yf, shape = stream)) +
geom_point() 
ggplot(met, aes(y = lER, x = Tanom.iktCg, color = Yf)) +
geom_point() +
facet_grid(stream ~.)
ggplot(met, aes(y = lER, x = invKT.C.StMean, color = Yf, shape = stream)) +
geom_point(position = "jitter") 
ggplot(met, aes(y = lER, x = invKT.C.GMean, color = Yf, shape = stream)) +
geom_point(position = "jitter") 
ER models
Likely maximal model
g.er.PutGLOB0 <- gamm(lER ~ treatment*invKT.C.StMean + s(meanQds.lcg, by = stream) , REML = FALSE, data = met)
g.er.PutGLOB1 <- gamm(lER ~ treatment*invKT.C.StMean + s(meanQds.lcg, by = stream) +
s(lGPP.cg, by = stream), REML = FALSE, data = met)
g.er.PutGLOB1a <- gamm(lER ~ treatment*invKT.C.StMean + s(meanQds.lcg, by = stream), REML = FALSE, data = met)
# this is perhaps the most intuitive, but it doesn't work
# g.er.PutGLOB1 <- gamm(lER ~ s(invKT.C.StMean, by = treatment) + s(meanQds.lc, by = stream) +
# s(LightPerDay.lc, by = stream) + s(Tanom.iktC, by = stream), REML = FALSE, data = met)model.sel(g.er.PutGLOB0$lme, g.er.PutGLOB1$lme, g.er.PutGLOB1a$lme)## Model selection table
## X family
## g.er.PutGLOB1$lme + gaussian(identity)
## g.er.PutGLOB0$lme + gaussian(identity)
## g.er.PutGLOB1a$lme + gaussian(identity)
## random
## g.er.PutGLOB1$lme X-g+X.0-g.0%i%g+X.1-g.1%i%g.0%i%g+X.2-g.2%i%g.1%i%g.0%i%g+X.3-g.3%i%g.2%i%g.1%i%g.0%i%g+X.4-g.4%i%g.3%i%g.2%i%g.1%i%g.0%i%g+X.5-g.5%i%g.4%i%g.3%i%g.2%i%g.1%i%g.0%i%g+X.6-g.6%i%g.5%i%g.4%i%g.3%i%g.2%i%g.1%i%g.0%i%g
## g.er.PutGLOB0$lme X-g+X.0-g.0%i%g+X.1-g.1%i%g.0%i%g+X.2-g.2%i%g.1%i%g.0%i%g
## g.er.PutGLOB1a$lme X-g+X.0-g.0%i%g+X.1-g.1%i%g.0%i%g+X.2-g.2%i%g.1%i%g.0%i%g
## df logLik AICc delta weight
## g.er.PutGLOB1$lme 23 -415.529 878.7 0.00 1
## g.er.PutGLOB0$lme 15 -610.288 1251.3 372.55 0
## g.er.PutGLOB1a$lme 15 -610.288 1251.3 372.55 0
## Models ranked by AICc(x)
## Random terms:
## X-g = 'Xr - 1 | g'
## X.0-g.0%i%g = 'Xr.0 - 1 | g.0 %in% g'
## X.1-g.1%i%g.0%i%g = 'Xr.1 - 1 | g.1 %in% g.0 %in% g'
## X.2-g.2%i%g.1%i%g.0%i%g = 'Xr.2 - 1 | g.2 %in% g.1 %in% g.0 %in% g'
## X.3-g.3%i%g.2%i%g.1%i%g.0%i%g = 'Xr.3 - 1 | g.3 %in% g.2 %in% g.1 %in% g.0 %in% g'
## X.4-g.4%i%g.3%i%g.2%i%g.1%i%g.0%i%g = 'Xr.4 - 1 | g.4 %in% g.3 %in% g.2 %in% g.1 %in% g.0 %in% g'
## X.5-g.5%i%g.4%i%g.3%i%g.2%i%g.1%i%g.0%i%g = 'Xr.5 - 1 | g.5 %in% g.4 %in% g.3 %in% g.2 %in% g.1 %in% g.0 %in% g'
## X.6-g.6%i%g.5%i%g.4%i%g.3%i%g.2%i%g.1%i%g.0%i%g = 'Xr.6 - 1 | g.6 %in% g.5 %in% g.4 %in% g.3 %in% g.2 %in% g.1 %in% g.0 %in% g'summary(g.er.PutGLOB1$gam)##
## Family: gaussian
## Link function: identity
##
## Formula:
## lER ~ treatment * invKT.C.StMean + s(meanQds.lcg, by = stream) +
## s(lGPP.cg, by = stream)
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.74379 0.07204 10.325 < 2e-16 ***
## treatmentnitrogen 1.69449 0.06688 25.337 < 2e-16 ***
## treatmentphosphorus 1.49621 0.05325 28.095 < 2e-16 ***
## invKT.C.StMean -0.57992 0.14679 -3.951 8.64e-05 ***
## treatmentnitrogen:invKT.C.StMean 0.41608 0.15170 2.743 0.00626 **
## treatmentphosphorus:invKT.C.StMean 0.57231 0.12917 4.431 1.10e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(meanQds.lcg):streamst11U 1.000 1.000 0.953 0.329
## s(meanQds.lcg):streamst18 5.023 5.023 12.240 1.96e-11 ***
## s(meanQds.lcg):streamst6 1.000 1.000 17.223 3.75e-05 ***
## s(meanQds.lcg):streamst9 3.214 3.214 22.691 1.71e-14 ***
## s(lGPP.cg):streamst11U 1.000 1.000 0.934 0.334
## s(lGPP.cg):streamst18 4.015 4.015 73.572 < 2e-16 ***
## s(lGPP.cg):streamst6 1.000 1.000 123.824 < 2e-16 ***
## s(lGPP.cg):streamst9 2.710 2.710 47.278 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.853
## Scale est. = 0.18624 n = 677g.er.PutGLOB0.b <- getViz(g.er.PutGLOB1$gam)check(g.er.PutGLOB0.b)## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
So autocorrelation
# https://www.kaggle.com/timib1203/managing-many-models-with-tidyr-purrr-and-broom
met.PutGLOB <- met %>%
mutate(g.erRes = resid(g.er.PutGLOB1$lme, type = "normalized")) %>%
select(streamTreat, g.erRes) %>%
group_by(streamTreat) %>%
nest() %>%
mutate(acf = map(.x = data, .f = function(d) ggAcf(d$g.erRes) + labs(title = streamTreat)))grid.arrange(grobs = met.PutGLOB$acf, ncol = 3)
Evaluate different corr structures
g.er.AC_1.0 <- gamm(lER ~ treatment*invKT.C.StMean + s(meanQds.lcg, by = stream),
correlation = corARMA(p = 1,q = 0, form = ~ doy | streamTreat), REML = FALSE, data = met)
g.er.AC_1.1 <- gamm(lER ~ treatment*invKT.C.StMean + s(meanQds.lcg, by = stream),
correlation = corARMA(p = 1,q = 1, form = ~ doy | streamTreat), REML = FALSE, data = met)
g.er.AC_1.2 <- gamm(lER ~ treatment*invKT.C.StMean + s(meanQds.lcg, by = stream),
correlation = corARMA(p = 1,q = 2, form = ~ doy | streamTreat), REML = FALSE, data = met)
g.er.AC_2.0 <- gamm(lER ~ treatment*invKT.C.StMean + s(meanQds.lcg, by = stream),
correlation = corARMA(p = 2,q = 0, form = ~ doy | streamTreat), REML = FALSE, data = met)
g.er.AC_2.1 <- gamm(lER ~ treatment*invKT.C.StMean + s(meanQds.lcg, by = stream),
correlation = corARMA(p = 2,q = 1, form = ~ doy | streamTreat), REML = FALSE, data = met)
g.er.AC_2.2 <- gamm(lER ~ treatment*invKT.C.StMean + s(meanQds.lcg, by = stream),
correlation = corARMA(p = 2,q = 2, form = ~ doy | streamTreat), REML = FALSE, data = met)
# doesn't work
g.er.AC_3.0 <- gamm(lER ~ treatment*invKT.C.StMean + s(meanQds.lcg, by = stream),
correlation = corARMA(p = 3,q = 0, form = ~ doy | streamTreat), REML = FALSE, data = met)
g.er.AC_3.1 <- gamm(lER ~ treatment*invKT.C.StMean + s(meanQds.lcg, by = stream),
correlation = corARMA(p = 3,q = 1, form = ~ doy | streamTreat), REML = FALSE, data = met)
#
g.er.AC_3.2 <- gamm(lER ~ treatment*invKT.C.StMean + s(meanQds.lcg, by = stream),
correlation = corARMA(p = 3,q = 2, form = ~ doy | streamTreat), REML = FALSE, data = met)Model selection
model.sel(g.er.AC_1.0$lme, g.er.AC_1.1$lme, g.er.AC_1.2$lme, g.er.AC_2.0$lme,
g.er.AC_2.1$lme, g.er.AC_2.2$lme, g.er.AC_3.0$lme, g.er.AC_3.1$lme, g.er.AC_3.2$lme) ## Model selection table
## X family df logLik AICc delta weight
## g.er.AC_3.2$lme + gaussian(identity) 20 -35.365 112.0 0.00 0.933
## g.er.AC_3.0$lme + gaussian(identity) 18 -41.171 119.4 7.37 0.023
## g.er.AC_3.1$lme + gaussian(identity) 19 -40.472 120.1 8.09 0.016
## g.er.AC_1.2$lme + gaussian(identity) 18 -41.761 120.6 8.55 0.013
## g.er.AC_2.2$lme + gaussian(identity) 19 -40.709 120.6 8.56 0.013
## g.er.AC_1.1$lme + gaussian(identity) 17 -45.297 125.5 13.51 0.001
## g.er.AC_2.0$lme + gaussian(identity) 17 -46.246 127.4 15.41 0.000
## g.er.AC_1.0$lme + gaussian(identity) 16 -48.863 130.5 18.54 0.000
## g.er.AC_2.1$lme + gaussian(identity) 18 -48.438 133.9 21.91 0.000
## Models ranked by AICc(x)
## Random terms (all models):
## 'Xr - 1 | g', 'Xr.0 - 1 | g.0 %in% g', 'Xr.1 - 1 | g.1 %in% g.0 %in% g', 'Xr.2 - 1 | g.2 %in% g.1 %in% g.0 %in% g'Much better. S6p, s6n, s9a, s9p, s11u still have some minor problems
# https://www.kaggle.com/timib1203/managing-many-models-with-tidyr-purrr-and-broom
# https://jroy042.github.io/nonlinear/week4.html
# https://stackoverflow.com/questions/47586110/autocorrelation-in-generalized-additive-models-gam
g.er.AC_best.acf <- met %>%
mutate(g.erRes = resid(g.er.AC_3.2$lme, type = "normalized")) %>%
select(streamTreat, g.erRes) %>%
group_by(streamTreat) %>%
nest() %>%
mutate(acf = map(.x = data, .f = function(d) ggAcf(d$g.erRes) + labs(title = streamTreat)))grid.arrange(grobs = g.er.AC_best.acf$acf, ncol = 3)
Let’s do some model selection. Thinking about this as our base model has Q, temp, and light. I’m not going to obsess about selecting those. Instead, I’ll just focus on Annual Temp and treatment
going to use a simpler one, problems fitting models with 3.0
BestERCor <- corARMA(p = 3,q = 2, form = ~ doy | streamTreat)# did quick check to see if I need T anom by stream - yes massive improvement in AIC
# also did check to see if I needed GPP x stream - nope!
# GLOBAL CENTERED Q
g.er.NutXT <- gamm(lER ~ treatment*invKT.C.StMean + s(meanQds.lcg, by = stream),
correlation = BestERCor, REML = FALSE, data = met)
g.er.NutpT <- gamm(lER ~ treatment + invKT.C.StMean + s(meanQds.lcg, by = stream),
correlation = BestERCor, REML = FALSE, data = met)
g.er.NutpTs <- gamm(lER ~ treatment + invKT.C.StMean + s(meanQds.lcs, by = stream),
correlation = BestERCor, REML = FALSE, data = met)
g.er.NutpTb <- gamm(lER ~ treatment + invKT.C.StMean + s(meanQds.lcg),
correlation = BestERCor, REML = FALSE, data = met)
g.er.NutXst <- gamm(lER ~ treatment*stream + s(meanQds.lcg, by = stream),
correlation = BestERCor, REML = FALSE, data = met)
g.er.Nutpst <- gamm(lER ~ treatment + stream + s(meanQds.lcg, by = stream),
correlation = BestERCor, REML = FALSE, data = met)
g.er.Nut <- gamm(lER ~ treatment + s(meanQds.lcg, by = stream),
correlation = BestERCor, REML = FALSE, data = met)
g.er.T <- gamm(lER ~ invKT.C.StMean + s(meanQds.lcg, by = stream),
correlation = BestERCor, REML = FALSE, data = met)
g.er.st <- gamm(lER ~ stream + s(meanQds.lcg, by = stream),
correlation = BestERCor, REML = FALSE, data = met)
g.er.null <- gamm(lER ~ s(meanQds.lcg, by = stream),
correlation = BestERCor, REML = FALSE, data = met)
# STREAM CENTERED Q
g.er.NutXTs <- gamm(lER ~ treatment*invKT.C.StMean + s(meanQds.lcs, by = stream),
correlation = BestERCor, REML = FALSE, data = met)
g.er.NutpTs <- gamm(lER ~ treatment + invKT.C.StMean + s(meanQds.lcs, by = stream),
correlation = BestERCor, REML = FALSE, data = met)
g.er.NutXsts <- gamm(lER ~ treatment*stream + s(meanQds.lcs, by = stream),
correlation = BestERCor, REML = FALSE, data = met)
g.er.Nutpsts <- gamm(lER ~ treatment + stream + s(meanQds.lcs, by = stream),
correlation = BestERCor, REML = FALSE, data = met)
g.er.Nuts <- gamm(lER ~ treatment + s(meanQds.lcs, by = stream),
correlation = BestERCor, REML = FALSE, data = met)
g.er.Ts <- gamm(lER ~ invKT.C.StMean + s(meanQds.lcs, by = stream),
correlation = BestERCor, REML = FALSE, data = met)
g.er.sts <- gamm(lER ~ stream + s(meanQds.lcs, by = stream),
correlation = BestERCor, REML = FALSE, data = met)
g.er.nulls <- gamm(lER ~ s(meanQds.lcs, by = stream),
correlation = BestERCor, REML = FALSE, data = met) Probably a tie between treatment and treatment + temp But the model with a stream centered Q is wicked better
model.sel(g.er.NutXT$lme, g.er.NutpT$lme, g.er.NutXst$lme, g.er.Nutpst$lme, g.er.Nut$lme, g.er.T$lme, g.er.st$lme, g.er.null$lme, g.er.NutpTb$lme,
g.er.NutXTs$lme, g.er.NutpTs$lme, g.er.NutXsts$lme, g.er.Nutpsts$lme, g.er.Nuts$lme, g.er.Ts$lme, g.er.sts$lme, g.er.nulls$lme)## Model selection table
## X family
## g.er.NutpTs$lme + gaussian(identity)
## g.er.Nutpsts$lme + gaussian(identity)
## g.er.NutXsts$lme + gaussian(identity)
## g.er.NutXTs$lme + gaussian(identity)
## g.er.Nuts$lme + gaussian(identity)
## g.er.Nut$lme + gaussian(identity)
## g.er.NutpT$lme + gaussian(identity)
## g.er.Nutpst$lme + gaussian(identity)
## g.er.NutXT$lme + gaussian(identity)
## g.er.NutXst$lme + gaussian(identity)
## g.er.nulls$lme + gaussian(identity)
## g.er.sts$lme + gaussian(identity)
## g.er.null$lme + gaussian(identity)
## g.er.T$lme + gaussian(identity)
## g.er.Ts$lme + gaussian(identity)
## g.er.st$lme + gaussian(identity)
## g.er.NutpTb$lme + gaussian(identity)
## random df
## g.er.NutpTs$lme X-g+X.0-g.0%i%g+X.1-g.1%i%g.0%i%g+X.2-g.2%i%g.1%i%g.0%i%g 18
## g.er.Nutpsts$lme X-g+X.0-g.0%i%g+X.1-g.1%i%g.0%i%g+X.2-g.2%i%g.1%i%g.0%i%g 20
## g.er.NutXsts$lme X-g+X.0-g.0%i%g+X.1-g.1%i%g.0%i%g+X.2-g.2%i%g.1%i%g.0%i%g 26
## g.er.NutXTs$lme X-g+X.0-g.0%i%g+X.1-g.1%i%g.0%i%g+X.2-g.2%i%g.1%i%g.0%i%g 20
## g.er.Nuts$lme X-g+X.0-g.0%i%g+X.1-g.1%i%g.0%i%g+X.2-g.2%i%g.1%i%g.0%i%g 17
## g.er.Nut$lme X-g+X.0-g.0%i%g+X.1-g.1%i%g.0%i%g+X.2-g.2%i%g.1%i%g.0%i%g 17
## g.er.NutpT$lme X-g+X.0-g.0%i%g+X.1-g.1%i%g.0%i%g+X.2-g.2%i%g.1%i%g.0%i%g 18
## g.er.Nutpst$lme X-g+X.0-g.0%i%g+X.1-g.1%i%g.0%i%g+X.2-g.2%i%g.1%i%g.0%i%g 20
## g.er.NutXT$lme X-g+X.0-g.0%i%g+X.1-g.1%i%g.0%i%g+X.2-g.2%i%g.1%i%g.0%i%g 20
## g.er.NutXst$lme X-g+X.0-g.0%i%g+X.1-g.1%i%g.0%i%g+X.2-g.2%i%g.1%i%g.0%i%g 26
## g.er.nulls$lme X-g+X.0-g.0%i%g+X.1-g.1%i%g.0%i%g+X.2-g.2%i%g.1%i%g.0%i%g 15
## g.er.sts$lme X-g+X.0-g.0%i%g+X.1-g.1%i%g.0%i%g+X.2-g.2%i%g.1%i%g.0%i%g 18
## g.er.null$lme X-g+X.0-g.0%i%g+X.1-g.1%i%g.0%i%g+X.2-g.2%i%g.1%i%g.0%i%g 15
## g.er.T$lme X-g+X.0-g.0%i%g+X.1-g.1%i%g.0%i%g+X.2-g.2%i%g.1%i%g.0%i%g 16
## g.er.Ts$lme X-g+X.0-g.0%i%g+X.1-g.1%i%g.0%i%g+X.2-g.2%i%g.1%i%g.0%i%g 16
## g.er.st$lme X-g+X.0-g.0%i%g+X.1-g.1%i%g.0%i%g+X.2-g.2%i%g.1%i%g.0%i%g 18
## g.er.NutpTb$lme X-g 12
## logLik AICc delta weight
## g.er.NutpTs$lme -32.672 102.4 0.00 0.447
## g.er.Nutpsts$lme -31.127 103.5 1.15 0.251
## g.er.NutXsts$lme -25.728 105.6 3.23 0.089
## g.er.NutXTs$lme -32.594 106.5 4.08 0.058
## g.er.Nuts$lme -35.776 106.5 4.10 0.058
## g.er.Nut$lme -35.840 106.6 4.22 0.054
## g.er.NutpT$lme -35.463 108.0 5.58 0.027
## g.er.Nutpst$lme -34.297 109.9 7.49 0.011
## g.er.NutXT$lme -35.365 112.0 9.63 0.004
## g.er.NutXst$lme -29.408 113.0 10.59 0.002
## g.er.nulls$lme -43.347 117.4 15.04 0.000
## g.er.sts$lme -42.082 121.2 18.82 0.000
## g.er.null$lme -45.494 121.7 19.33 0.000
## g.er.T$lme -45.321 123.5 21.08 0.000
## g.er.Ts$lme -46.578 126.0 23.60 0.000
## g.er.st$lme -44.544 126.1 23.74 0.000
## g.er.NutpTb$lme -59.847 144.2 41.78 0.000
## Models ranked by AICc(x)
## Random terms:
## X-g = 'Xr - 1 | g'
## X.0-g.0%i%g = 'Xr.0 - 1 | g.0 %in% g'
## X.1-g.1%i%g.0%i%g = 'Xr.1 - 1 | g.1 %in% g.0 %in% g'
## X.2-g.2%i%g.1%i%g.0%i%g = 'Xr.2 - 1 | g.2 %in% g.1 %in% g.0 %in% g'summary(g.er.NutpTs$gam)##
## Family: gaussian
## Link function: identity
##
## Formula:
## lER ~ treatment + invKT.C.StMean + s(meanQds.lcs, by = stream)
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.8139 0.2218 3.669 0.000263 ***
## treatmentnitrogen 1.7862 0.3011 5.933 4.81e-09 ***
## treatmentphosphorus 1.3471 0.2994 4.500 8.05e-06 ***
## invKT.C.StMean -1.0141 0.3709 -2.734 0.006426 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(meanQds.lcs):streamst11U 2.484 2.484 7.461 0.000658 ***
## s(meanQds.lcs):streamst18 3.597 3.597 4.288 0.002139 **
## s(meanQds.lcs):streamst6 4.039 4.039 11.987 1.71e-09 ***
## s(meanQds.lcs):streamst9 5.030 5.030 20.046 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.609
## Scale est. = 0.46633 n = 677g.er.NutpT.b <- getViz(g.er.NutpTs$gam)This model violates some shit!!!
check(g.er.NutpT.b)## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
hist(resid(g.er.NutpT$lme, type = "normalized"))
g.er.NutpT.acf <- met %>%
mutate(g.erRes = resid(g.er.NutpTs$lme, type = "normalized")) %>%
select(streamTreat, g.erRes) %>%
group_by(streamTreat) %>%
nest() %>%
mutate(acf = map(.x = data, .f = function(d) ggAcf(d$g.erRes) + labs(title = streamTreat)))grid.arrange(grobs = g.er.NutpT.acf$acf, ncol = 3)
met$BestMod <- predict.gam(g.er.NutpTs$gam, met, type = "response")ggplot(met, aes(y = lER, x = BestMod, color = treatment)) +
geom_point()+
facet_grid(stream ~.) +
geom_abline(intercept = 0, slope =1)
print(plot(g.er.NutpT.b, allTerms = T), pages = 1)
In the splines, this shows all the data!
blah <- plot(g.er.NutpT.b, allTerms = 1) +
l_points(color = "blue") + l_fitLine(linetype = 3) + l_fitContour() +
l_ciLine(colour = 2) + l_ciBar() + l_fitPoints(size = 1, col = 2) + theme_get() + labs(title = NULL)
print(blah, pages = 1)
save.image("01_ERanalysis_rdat")
# load("01_ERanalysis_rdat")